%% This file generates results for the phrase "Among farmers whose expected
%  number of trees is above the reward threshold at time of take up, XX% do
%  not reach the threshold bc their ex post number of trees is lower than
%  their ex ante number, once new information is accounted for 

% See exAnteexPostTrees_compare.xls

clear
load resultsFMIN.mat

clearvars snPass


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% 1. Getting ready for graphs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

% Cost parameters (in thousand Kwacha);
cp.alpha  = 48;           % Public benefit per tree;
cp.eta    = 1.2;          % Cost of public funds;
cp.mon    = 35;           % Monitoring cost;
cp.cb     = (1/(1+0.03)); % Cost of borrowing;

% Run parameters, tolerance is unchanged
rp.parl    = 0;
rp.m       = 200;
rp.k       = 1;
rp.lam     = 0;
rp.theseed = rp.theseed+1000;


rng(rp.theseed);
snPass.s       = rand(obs*(3*rp.k + rp.m),1);

% These two lines create a (k x m x 51) matrix, where the first floor is a
% (k x m) matrix of ones, the second a (k x m) matrix of twos, up to 51;
% Nv is just N repeated k x 51 times;
% the fourth and fifth are true false matrices for N vs. Nbar;
Ntr            = permute(N,[3,1,2]);
snPass.N3D     = Ntr(ones(rp.k,1),ones(rp.m,1),:);
snPass.Nv      = N(ones(rp.k,1),:);
snPass.Ngeq3D  = (snPass.N3D>=Nbar);
snPass.N03D    = (snPass.N3D>0);
snPass.Ngeqv   = (snPass.Nv>=Nbar);
snPass.N0v     = (snPass.Nv>0);


% Must reset several values, setting to 0: F_surp (F_treat), F_moni, and
% F_shift -- This removes this context specific shifters. Further, must
% redefine Nmatrices because of new k/m dimensions

theta_hat0 = coefftran([theta_hat2(1:6),0,0,0],0); % Ensures F_* parameters = 0
theta_hat1 = coefftran([theta_hat2(1:4),-50,theta_hat2(6),0,0,0],0); % Ensures F_* parameters = 0


TG0   = zeros(obs,1);      % Alternate treatment vector, set = 1
MON0  = zeros(obs,1);      % Alternate monitoring vector, set = 0


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2a) Tree Survival and comparisons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[PartProb,Npdf,Nobs,Prof1,Prof0,EN,OV,Prof0_C] = Full_Simulate(theta_hat0,DP,A,R,TG0,MON0,snPass,rp);


%% Relevent
1314-sum(PartProb) % Particpation
EN = round(EN);

%%

%Row1
    [sum(PartProb==1 & EN==0 & Nobs==0),...
    sum(PartProb==1 & EN==0 & Nobs>0 & Nobs<35),...
    sum(PartProb==1 & EN==0 & Nobs==35),...
    sum(PartProb==1 & EN==0 & Nobs>35 & Nobs<50),...
    sum(PartProb==1 & EN==0 & Nobs==50);...
 %Row2
    sum(PartProb==1 & EN>0 & EN<35 & Nobs==0),...
    sum(PartProb==1 & EN>0 & EN<35 & Nobs>0 & Nobs<35),...
    sum(PartProb==1 & EN>0 & EN<35 & Nobs==35),...
    sum(PartProb==1 & EN>0 & EN<35 & Nobs>35 & Nobs<50),...
    sum(PartProb==1 & EN>0 & EN<35 & Nobs==50);...   
 %Row3
    sum(PartProb==1 & EN==35 & Nobs==0),...
    sum(PartProb==1 & EN==35 & Nobs>0 & Nobs<35),...
    sum(PartProb==1 & EN==35 & Nobs==35),...
    sum(PartProb==1 & EN==35 & Nobs>35 & Nobs<50),...
    sum(PartProb==1 & EN==35 & Nobs==50);...
 %Row4
    sum(PartProb==1 & EN>35 & EN<50 & Nobs==0),...
    sum(PartProb==1 & EN>35 & EN<50 & Nobs>0 & Nobs<35),...
    sum(PartProb==1 & EN>35 & EN<50 & Nobs==35),...
    sum(PartProb==1 & EN>35 & EN<50 & Nobs>35 & Nobs<50),...
    sum(PartProb==1 & EN>35 & EN<50 & Nobs==50);...
 %Row5
    sum(PartProb==1 & EN==50 & Nobs==0),...
    sum(PartProb==1 & EN==50 & Nobs>0 & Nobs<35),...
    sum(PartProb==1 & EN==50 & Nobs==35),...
    sum(PartProb==1 & EN==50 & Nobs>35 & Nobs<50),...
    sum(PartProb==1 & EN==50 & Nobs==50)]


%%

sum(PartProb==1 & Nobs>EN)
sum(PartProb==1 & Nobs==EN)
sum(PartProb==1 & Nobs<EN)

sum(PartProb==1 & Nobs>0 & Nobs>EN)
sum(PartProb==1 & Nobs>0 & Nobs==EN)
sum(PartProb==1 & Nobs>0 & Nobs<EN)
%%
sum(PartProb==1 & EN>0 & Nobs>EN)

%%
%sum(PartProb==1 & Nobs>0 & Nobs<35 & EN==35) / sum(PartProb==1 & Nobs>0 & EN==35)
%sum(PartProb==1 & Nobs>0 & Nobs<35 & EN>=35) / sum(PartProb==1 & Nobs>0 & EN>=35)
%sum(PartProb==1 & Nobs==0 & EN==35)
%sum(PartProb==1 & Nobs==0 & EN>=35)

sum(PartProb==1 & Nobs<35 & EN==35) / sum(PartProb==1 & Nobs>0 & EN==35)
sum(PartProb==1 & Nobs<35 & EN>=35) / sum(PartProb==1 & Nobs>0 & EN>=35)

